#include "cubaturs.h"

#define EM_E0 GSL_CONST_MKSA_VACUUM_PERMITTIVITY //8.854187817e-12
#define EM_M0 GSL_CONST_MKSA_VACUUM_PERMEABILITY //1.25663706144e-6	//electric&magnetic constants
	
	int j;
	
void Cubaturs::calc(int area, int method)
{
	for(int j=0;j<9; j++)	//Инициализация
	{
		ss[j]=0I;	//Инициализация
		RES[j]=ss[j];
	}
        double ee=1,mm=1,d;
	__complex__ double di;
	Integral=0I;
	switch (area)
	{
		case 0:
			ro0=0; fis0=0;
                        dq=a/N; dp=2*M_PI/N;
                        ee=e1*EM_E0; mm=m1*EM_M0;
			break;
		case 1 :
			ro0=0; fis0=0;
                        dq=a/N; dp=2*M_PI/N;
                        ee=e1*EM_E0; mm=m1*EM_M0;
			break;
		case 2 :
			ro0=a; fis0=0;
                        dq=(b-a)/N; dp=2*M_PI/N;
                        ee=e2*EM_E0; mm=m2*EM_M0;
			break;
		case 3 :
			break;
		case 4 :
			ro0=0; fis0=0;
                        dq=a/N; dp=2*M_PI/N;
                        ee=-1*m1*EM_M0; mm=-1*e1*EM_E0;
			break;
		case 5 :
			ro0=0; fis0=0;
                        dq=a/N; dp=2*M_PI/N;
                        ee=e1*EM_E0; mm=m1*EM_M0;
			break;
		case 6 :
			ro0=0; fis0=0;
                        dq=a/N; dp=2*M_PI/N;
                        ee=e1*EM_E0; mm=m1*EM_M0;
			break;
		case 7 :
			ro0=0; fis0=0;
                        dq=a/N; dp=2*M_PI/N;
                        ee=e1*EM_E0; mm=m1*EM_M0;
                        break;
                case 8 :
                        ro0=0; fis0=0;
                        dq=a/N; dp=2*M_PI/N;
                        ee=e1*EM_E0; mm=m1*EM_M0;
                        break;
	}
	bb[0]=0.2113248654; 
	bb[1]=0.7886751346;
	ii2=0I;
	ia=0I; ib=0I; ic=0I; id=0I;
	for(int p=0;p<N;p++)
	{
		for(int q=0;q<N;q++)
		{
			for(int s1=0;s1<2;s1++)
			{
				fis=fis0+bb[s1]*dp;
				for(int s2=0;s2<2;s2++)
				{
					ro=ro0+bb[s2]*dq;
					switch (area)
					{
					case 0 :
						Ero=cos(fis)+0I;
						Efi=-1*sin(fis)+0I;
						Hro=sqrt(ee/mm)*sin(fis)+0I;
						Hfi=sqrt(ee/mm)*cos(fis)+0I;
                                                break;
					case 1 :
						Ero=Ero1(ro,fis);
						Efi=Efi1(ro,fis);
						Hro=Hro1(ro,fis);
						Hfi=Hfi1(ro,fis);
  						break;
					case 2 :
						Ero=Ero2(ro,fis);
						Efi=Efi2(ro,fis);
						Hro=Hro2(ro,fis);
						Hfi=Hfi2(ro,fis);
						break;
					case 3 : //- reserved for main structure
						break;
					case 4 :	// Поле H для равноамплитудного возбуждения
						Hro=sin(fis)+0I;
						Hfi=cos(fis)+0I;
						Ero=sqrt(mm/ee)*cos(fis)+0I;
						Efi=-1*sqrt(mm/ee)*sin(fis)+0I;
						break;
					case 5 :	//радиальноуменьшающееся возбуждение
						Ero=cos(fis)*(1-pow(ro/a,2))+0I;
						Efi=-1*sin(fis)*(1-pow(ro/a,2))+0I;
						Hro=sqrt(ee/mm)*sin(fis)*(1-pow(ro/a,2))+0I;
						Hfi=sqrt(ee/mm)*cos(fis)*(1-pow(ro/a,2))+0I;
						break;
					case 6 :
                                                Ero=-1*1I*beta/alfa1*js1(ro);    //?
						Efi=0I;
						Hro=0I;
                                                Hfi=-1*1I*w*ee/alfa1*js1(ro);    //?
						break;
					case 7 :
						Ero=0I;
                                                Efi=1I*w*mm/alfa1*js1(ro);       //?
                                                Hro=-1*1I*beta/alfa1*js1(ro);    //?
						Hfi=0I;
                                                break;
                                        case 8 :
                                                Ero=1I*w*mm/pow(alfa1,2)/ro*gsl_sf_bessel_J1(alfa1*ro)*sin(fis);
                                                Efi=1I*w*mm/alfa1*js1(ro)*cos(fis);
                                                Hro=-1*1I*beta/alfa1*js1(ro)*cos(fis);
                                                Hfi=1I*beta/pow(alfa1,2)/ro*gsl_sf_bessel_J1(alfa1*ro)*sin(fis);
                                                break;
					}
					switch (method)
					{
					case 0 :
					//Численное интегрирование по строгому алгоритму
						rr=sqrt(pow(R,2)+pow(ro,2)-2*R*ro*cosg(fis));
                                                ri=(1/rr+1I*k);
					//ErR
						ii1[0]=Ero*ro/rr*ri;
						ii1[0]*=ro*cexp(rr);
						ss[0]+=ii1[0];
					//ErT
						ii1[1]=1/rr*ri*(R*(Ero*cos(fi-fis)+Efi*sin(fi-fis))-ro*Ero*sin(teta));
						ii1[1]*=ro*cexp(rr);
						ss[1]+=ii1[1];
					//ErF
						ii1[2]=ri/rr*(Ero*cosb(fis)-Efi*cosd(fis));
						ii1[2]*=ro*cexp(rr);
						ss[2]+=ii1[2];
					//Public multiplier for grad_div
						ic=3/rr*ri-pow(k,2);
						ic*=(R*(Hro*cosa(fis)-Hfi*cosg(fis))+ro*Hfi)/pow(rr,2);
					//EgR
						ia=(R-ro*cosg(fis))*ic;
						ib=ri/rr*(Hro*cosa(fis)-Hfi*cosg(fis));
						ii1[3]=ia-ib;
						ii1[3]*=ro*cexp(rr);
						ss[3]+=ii1[3];
					//EgT
						ia=ro*cosd(fis)*ic;
						ib=ri/rr*(Hro*cosb(fis)-Hfi*cosd(fis));
						ii1[4]=ia+ib;
						ii1[4]*=ro*cexp(rr);
						ss[4]+=ii1[4];
					//EgF
						ia=ro*sin(fi-fis)*ic;
						ib=ri/rr*(Hro*cos(fi-fis)+Hfi*sin(fi-fis));
						ii1[5]=ia-ib;
						ii1[5]*=ro*cexp(rr);
						ss[5]+=ii1[5];
					//EjR
						ii1[6]=Hro*cosa(fis)-Hfi*cosg(fis);
						ii1[6]*=ro*cexp(rr);
						ss[6]+=ii1[6];
					//EjT
						ii1[7]=Hro*cosb(fis)-Hfi*cosd(fis);
						ii1[7]*=ro*cexp(rr);
						ss[7]+=ii1[7];
					//EjF
						ii1[8]=Hro*cos(fi-fis)+Hfi*sin(fi-fis);
						ii1[8]*=ro*cexp(rr);
						ss[8]+=ii1[8];
						break;
					case 1 :
						rr=sqrt(pow(R,2)+pow(ro,2)-2*R*ro*cosg(fis));
                                                ri=(1/rr+1I*k);
					//Строгий алгоритм для дальней зоны	
						d=R-ro*cosg(fis);
                                                di=(1/d+1I*k);
					//ErR
						ii1[0]=Ero*ro;
						ii1[0]*=ro*cexp0(ro,fis);
						ss[0]+=ii1[0];
					//ErT
						ii1[1]=ro*(Ero*cosa(fis)+Efi*cosg(fis))*sin(fi-fis);
						ii1[1]*=ro*cexp0(ro,fis);
						ss[1]+=ii1[1];
					//ErF
						ii1[2]=ro*(Ero*cosa(fis)+Efi*cosg(fis))*cosd(fis);
						ii1[2]*=ro*cexp0(ro,fis);
						ss[2]+=ii1[2];
					//EgR
                                                ii1[3]=(Hro*cosa(fis)-Hfi*cosg(fis))*(2-1I*k*ro*cosg(fis)-1I*k*ro*Hfi);
                                                ss[3]+=ii1[3];
					//EgT
						ia=(Hro*cosb(fis)-Hfi*cosd(fis))*(2+pow(k*ro,2)*cosa(fis)*cosg(fis));
						ib=pow(k*ro,2)*Hfi*cosd(fis);
						ii1[4]=ia-ib;
						ii1[4]*=ro*cexp0(ro,fis);
						ss[4]+=ii1[4];
					//EgF
                                                ia=1I*k*ro*cosa(fis)*sin(teta)*Hro+pow(k*ro,2)*cosb(fis)*cos(teta)*Hfi;
						ib=(Hro*cos(fi-fis)+Hfi*sin(fi-fis))*(2-pow(k*ro,2)*cosa(fis)*sin(teta));
						ii1[5]=ib-ia;
						ii1[5]*=ro*cexp0(ro,fis);
						ss[5]+=ii1[5];
					//EjR
						ii1[6]=Hro*cosa(fis)-Hfi*cosg(fis);
						ii1[6]*=ro*cexp0(ro,fis);
						ss[6]+=ii1[6];
					//EjT
						ii1[7]=Hro*cosb(fis)-Hfi*cosd(fis);
						ii1[7]*=ro*cexp0(ro,fis);
						ss[7]+=ii1[7];
					//EjF
						ii1[8]=Hro*cos(fi-fis)+Hfi*sin(fi-fis);
						ii1[8]*=ro*cexp0(ro,fis);
						ss[8]+=ii1[8]; 
						break;
					case 2 :
						rr=sqrt(pow(R,2)+pow(ro,2)+pow(Z,2)-2*R*ro*cos(fi-fis));
                                                ri=(1/rr+1I*k);
					//Строгий алгоритм для цилиндрической системы координат
					//ErR
						ii1[0]=(Ero*cos(fi-fis)+Efi*sin(fi-fis))*ri/rr;
						ii1[0]*=ro*cexp(rr);
						ss[0]+=ii1[0];
					//ErF
						ii1[1]=(Ero*sin(fi-fis)-Efi*cos(fi-fis))*ri/rr;
						ii1[1]*=ro*cexp(rr);
						ss[1]+=ii1[1];
					//ErZ
						ii1[2]=(ro*Ero-R*(Ero*cos(fi-fis)+Efi*sin(fi-fis)))*ri/rr;
						ii1[2]*=ro*cexp(rr);
						ss[2]+=ii1[2];
					//Public multiplier for grad_div
                                                ic=1/pow(rr,2)+ri*(2/rr+1I*k);
						ic*=(R*(Hro*sin(fi-fis)-Hfi*cos(fi-fis))+ro*Hfi)/pow(rr,2);
					//EgR
						ii1[3]=(R-ro*cos(fi-fis))*ic-ri/rr*(Hro*sin(fi-fis)-Hfi*cos(fi-fis));
						ii1[3]*=ro*cexp(rr);
						ss[3]+=ii1[3];
					//EgT
						ii1[4]=ro*sin(fi-fis)*ic-ri/rr*(Hro*cos(fi-fis)+Hfi*sin(fi-fis));
						ii1[4]*=ro*cexp(rr);
						ss[4]+=ii1[4];
					//EgF
						ii1[5]=ic;
						ii1[5]*=ro*cexp(rr);
						ss[5]+=ii1[5];
					//EjR
						ii1[6]=Hro*sin(fi-fis)-Hfi*cos(fi-fis);
						ii1[6]*=ro*cexp(rr);
						ss[6]+=ii1[6];
					//EjT
						ii1[7]=Hro*cos(fi-fis)+Hfi*sin(fi-fis);
						ii1[7]*=ro*cexp(rr);
						ss[7]+=ii1[7];
					//EjF
                                                ii1[8]*=1;
						ii1[8]*=ro*cexp(rr);
						ss[8]+=ii1[8];
						break;
					}
				}
			}
                        for(int j=0; j<9; j++) RES[j]+=ss[j];	//Прибавление результата по очередной кубатуре
			for(int j=0;j<9; j++) ss[j]=0I;	//инициализация нового узла
			ro0+=dq;		//инициализация нового узла
		}
		fis0+=dp; //инициализация новой строки кубатур
		ro0=0;	//возвращение каретки
	}
	for(int j=0; j<9; j++) RES[j]*=dp*dq/4;	//ending of Integer
	switch (method)
	{
		case 0 :
	//Численное интегрирование по строгому алгоритму
        RES[0]*=(cos(teta)/4/M_PI);
        RES[1]/=(4*M_PI);
        RES[2]*=(-1*R/4/M_PI);
        RES[3]/=(4*M_PI*1I*w*ee);
        RES[4]/=(-4*M_PI*1I*w*ee);
        RES[5]/=(4*M_PI*1I*w*ee);
        RES[6]*=(-1)*1I*w*mm/4/M_PI;
        RES[7]*=(-1)*1I*w*mm/4/M_PI;
        RES[8]*=(-1)*1I*w*mm/4/M_PI;
		break;
		case 1 :
	//Строгий алгоритм для дальней зоны
        RES[0]*=cexp(R)*1I*k*cos(teta)/(4*M_PI*R);
        RES[1]*=cexp(R)*1I*k/(-4*M_PI*R);
        RES[2]*=cexp(R)*1I*k/(-4*M_PI*R);
        RES[3]/=(4*M_PI*R*1I*w*ee);
        RES[4]*=cexp(R)/(4*M_PI*R*R*1I*w*ee);
        RES[5]*=cexp(R)/(4*M_PI*R*R*1I*w*ee);
        RES[6]*=cexp(R)*(-1)*1I*w*mm/4/M_PI;
        RES[7]*=cexp(R)*(-1)*1I*w*mm/4/M_PI;
        RES[8]*=cexp(R)*(-1)*1I*w*mm/4/M_PI;
		break;
		case 2:
	//Строгий алгоритм для цилиндрической системы координат
        RES[0]*=Z/4/M_PI;
        RES[1]*=Z/(-4*M_PI);
        RES[2]/=(4*M_PI);
        RES[3]/=(4*M_PI*1I*w*ee);
        RES[4]/=(4*M_PI*1I*w*ee);
        RES[5]*=Z/(4*M_PI*1I*w*ee);
        RES[6]*=(-1)*1I*w*mm/4/M_PI;
        RES[7]*=(-1)*1I*w*mm/4/M_PI;
        RES[8]*=(-1)*1I*w*mm/4/M_PI;
		break;
	}
}
